
### This script computes all robustness checks except weighting
### and rerunning the models with the observed value approach (both separate scripts)

library(haven)
library(sjPlot)
library(sjlabelled)
library(tidyverse)
library(gridExtra)
library(MASS)

##### Supplementary Appendix H: Robustness checks for Hypothesis 1 ####

### Table H2

# with first definition of activism
mH1_robset1_1 <- glm(polact1 ~ openST + conscST + extraST + agreeST + neurST +
                       age + edu + female + wealth + moscow,
                     family = "binomial", data = data1)

mH1_robset1_2 <- glm(polact1 ~ openST + conscST + extraST + agreeST + neurST +
                       age + edu + female + wealth + moscow + putin_sup,
                     family = "binomial", data = data1)

mH1_robset1_3 <- glm(polact1 ~ openST + conscST + extraST + agreeST + neurST +
                       age + edu + female + wealth + moscow + putin_app,
                     family = "binomial", data = data1)

mH1_robset1_4 <- glm(polact1 ~ openST + conscST + extraST + agreeST + neurST +
                       age + edu + female + wealth + moscow + poldisc,
                     family = "binomial", data = data1)

mH1_robset1_5 <- glm(polact1 ~ openST + conscST + extraST + agreeST + neurST +
                       age + edu + female + wealth + moscow + poldisc + putin_sup,
                     family = "binomial", data = data1)

mH1_robset1_6 <- glm(polact1 ~ openST + conscST + extraST + agreeST + neurST +
                       age + edu + female + wealth + moscow + poldisc + putin_app,
                     family = "binomial", data = data1)

# with alternative definition of activism
mH1_robset2_1 <- glm(polact2 ~ openST + conscST + extraST + agreeST + neurST +
                age + edu + female + wealth + moscow,
              family = "binomial", data = data1)

mH1_robset2_2 <- glm(polact2 ~ openST + conscST + extraST + agreeST + neurST +
                 age + edu + female + wealth + moscow + putin_sup,
               family = "binomial", data = data1)

mH1_robset2_3 <- glm(polact2 ~ openST + conscST + extraST + agreeST + neurST +
                 age + edu + female + wealth + moscow + putin_app,
               family = "binomial", data = data1)

mH1_robset2_4 <- glm(polact2 ~ openST + conscST + extraST + agreeST + neurST +
                       age + edu + female + wealth + moscow + poldisc,
                     family = "binomial", data = data1)

mH1_robset2_5 <- glm(polact2 ~ openST + conscST + extraST + agreeST + neurST +
                 age + edu + female + wealth + moscow + poldisc + putin_sup,
               family = "binomial", data = data1)

mH1_robset2_6 <- glm(polact2 ~ openST + conscST + extraST + agreeST + neurST +
                 age + edu + female + wealth + moscow + poldisc + putin_app,
               family = "binomial", data = data1)


H1_models_rob <- list(mH1_robset1_1, mH1_robset1_2, mH1_robset1_3, 
                       mH1_robset1_4, mH1_robset1_5, mH1_robset1_6,
                       mH1_robset2_1, mH1_robset2_2, mH1_robset2_3, 
                       mH1_robset2_4, mH1_robset2_5, mH1_robset2_6)
H1_modelnames_rob <- c("1", "2", "3", "4", "5", "6",
                        "7", "8", "9", "10", "11", "12")

tab_model(H1_models_rob,
          show.se = TRUE, show.ci = FALSE, 
          p.style = "stars", collapse.se = TRUE, show.ci50 = TRUE, show.aic = TRUE,
          file = "Table_H2.html")


##### Supplementary Appendix H: Robustness checks for Hypothesis 2 #####

### Table H3 on restricted data set ###

# restriction variable
data2$restrict <- ifelse(data2$edu>=3 & data2$wealth>=3,1,0)

# models
mH2_robset1_1 <- glm(navact ~ openST + conscST + extraST + agreeST + neurST,
                     family = "binomial", data = data2[data2$restrict==1,])

mH2_robset1_2 <- glm(navact ~ openST + conscST + extraST + agreeST + neurST +
                       age + edu + female + wealth + moscow,
                     family = "binomial", data = data2[data2$restrict==1,])

mH2_robset1_3 <- polr(factor(navact_scale) ~ openST + conscST + extraST + agreeST + neurST,
                      method = "logistic", data = data2[data2$restrict==1,])

mH2_robset1_4 <- polr(factor(navact_scale) ~ openST + conscST + extraST + agreeST + neurST +
                        age + edu + female + wealth + moscow,
                      method = "logistic", data = data2[data2$restrict==1,])

tab_model(mH2_robset1_1, mH2_robset1_2, mH2_robset1_3, mH2_robset1_4,
          show.se = TRUE, show.ci = FALSE, 
          p.style = "stars", collapse.se = TRUE, show.ci50 = TRUE, show.aic = TRUE,
          file = "Table_H3.html")


### Table H4 on full data set ###

mH2_robset2_1 <- glm(navact ~ openST + conscST + extraST + agreeST + neurST,
              family = "binomial", data = data2)

mH2_robset2_2 <- glm(navact ~ openST + conscST + extraST + agreeST + neurST +
                age + edu + female + wealth + moscow + poldisc,
              family = "binomial", data = data2)

mH2_robset2_3 <- polr(factor(navact_scale) ~ openST + conscST + extraST + agreeST + neurST,
               method = "logistic", data = data2)

mH2_robset2_4 <- polr(factor(navact_scale) ~ openST + conscST + extraST + agreeST + neurST +
                 age + edu + female + wealth + moscow + poldisc,
               method = "logistic", data = data2)

tab_model(mH2_robset2_1, mH2_robset2_2, mH2_robset2_3, mH2_robset2_4,
          show.se = TRUE, show.ci = FALSE, 
          p.style = "stars", collapse.se = TRUE, show.ci50 = TRUE, show.aic = TRUE,
          file = "Table_H4.html")


##### Supplementary Appendix H: Robustness checks for Hypothesis 3 within 2013 data #####

### create alternative subset variables to produce groups for case-control designs

data1$subset2 <- ifelse(data1$putin_app==0 & data1$polact1==1,1,
                                    ifelse(data1$putin_app==1 & data1$polact1==0,0,NA))

data1$subset3 <- ifelse(data1$putin_sup==0 & data1$polact2==1,1,
                        ifelse(data1$putin_sup==1 & data1$polact2==0,0,NA))

data1$subset4 <- ifelse(data1$putin_app==0 & data1$polact2==1,1,
                        ifelse(data1$putin_app==1 & data1$polact2==0,0,NA))


### Table H5: main models with added control variables (same case-control as main text)

mH3_robset1_1 <- glm(subset ~ openST + conscST + extraST + agreeST + neurST,
                     family = "binomial", data = data1)

mH3_robset1_2 <- glm(subset ~ openST + conscST + extraST + agreeST + neurST +
                       age + edu + female + wealth + moscow,
                     family = "binomial", data = data1)

mH3_robset1_3 <- glm(subset ~ openST + conscST + extraST + agreeST + neurST +
                       age + edu + female + wealth + moscow + poldiscST + econ + private,
                     family = "binomial", data = data1)

# models with alternative Putin support variable
mH3_robset2_1 <- glm(subset2 ~ openST + conscST + extraST + agreeST + neurST,
                     family = "binomial", data = data1)

mH3_robset2_2 <- glm(subset2 ~ openST + conscST + extraST + agreeST + neurST +
                       age + edu + female + wealth + moscow,
                     family = "binomial", data = data1)

mH3_robset2_3 <- glm(subset2 ~ openST + conscST + extraST + agreeST + neurST +
                       age + edu + female + wealth + moscow + poldiscST + econ + private,
                     family = "binomial", data = data1)

# models with alternative activism variable
mH3_robset3_1 <- glm(subset3 ~ openST + conscST + extraST + agreeST + neurST,
                     family = "binomial", data = data1)

mH3_robset3_2 <- glm(subset3 ~ openST + conscST + extraST + agreeST + neurST +
                       age + edu + female + wealth + moscow,
                     family = "binomial", data = data1)

mH3_robset3_3 <- glm(subset3 ~ openST + conscST + extraST + agreeST + neurST +
                       age + edu + female + wealth + moscow + poldiscST + econ + private,
                     family = "binomial", data = data1)

# models with alternative activism variable and alternative Putin support
mH3_robset4_1 <- glm(subset4 ~ openST + conscST + extraST + agreeST + neurST,
                 family = "binomial", data = data1)

mH3_robset4_2 <- glm(subset4 ~ openST + conscST + extraST + agreeST + neurST +
                   age + edu + female + wealth + moscow,
                 family = "binomial", data = data1)

mH3_robset4_3 <- glm(subset4 ~ openST + conscST + extraST + agreeST + neurST +
                   age + edu + female + wealth + moscow + poldiscST + econ + private,
                 family = "binomial", data = data1)

tab_model(mH3_robset1_1, mH3_robset1_2, mH3_robset1_3, 
          mH3_robset2_1, mH3_robset2_2, mH3_robset2_3,
          mH3_robset3_1, mH3_robset3_2, mH3_robset3_3,
          mH3_robset4_1, mH3_robset4_2, mH3_robset4_3,
          show.se = TRUE, show.ci = FALSE, 
          p.style = "stars", collapse.se = TRUE, show.ci50 = TRUE, show.aic = TRUE,
          file = "Table_H5.html")

##### Supplementary Appendix H: Robustness checks for Hypothesis 3 with pooled data #####

### create alternative subset variables to produce groups for case-control designs

pooled$subset2 <- ifelse(pooled$navact==1,1,
                         ifelse(pooled$putin_app==1 & pooled$polact1==0,0,NA))

pooled$subset3 <- ifelse(pooled$navact==1,1,
                         ifelse(pooled$putin_sup==1 & pooled$polact2==0,0,NA))

pooled$subset4 <- ifelse(pooled$navact==1,1,
                         ifelse(pooled$putin_app==1 & pooled$polact2==0,0,NA))


### Table H6: main models with added control variables (same case-control as main text)
mH3_robset5_1 <- glm(subset ~ openST + conscST + extraST + agreeST + neurST,
                     family = "binomial", data = pooled)

mH3_robset5_2 <- glm(subset ~ openST + conscST + extraST + agreeST + neurST +
                       age + edu + female + wealth + moscow,
                     family = "binomial", data = pooled)

mH3_robset5_3 <- glm(subset ~ openST + conscST + extraST + agreeST + neurST +
                       age + edu + female + wealth + moscow + poldiscST,
                     family = "binomial", data = pooled)

# models with alternative Putin support variable
mH3_robset6_1 <- glm(subset2 ~ openST + conscST + extraST + agreeST + neurST,
                     family = "binomial", data = pooled)

mH3_robset6_2 <- glm(subset2 ~ openST + conscST + extraST + agreeST + neurST +
                       age + edu + female + wealth + moscow,
                     family = "binomial", data = pooled)

mH3_robset6_3 <- glm(subset2 ~ openST + conscST + extraST + agreeST + neurST +
                       age + edu + female + wealth + moscow + poldiscST,
                     family = "binomial", data = pooled)

# models with alternative activism variable
mH3_robset7_1 <- glm(subset3 ~ openST + conscST + extraST + agreeST + neurST,
                     family = "binomial", data = pooled)

mH3_robset7_2 <- glm(subset3 ~ openST + conscST + extraST + agreeST + neurST +
                       age + edu + female + wealth + moscow,
                     family = "binomial", data = pooled)

mH3_robset7_3 <- glm(subset3 ~ openST + conscST + extraST + agreeST + neurST +
                       age + edu + female + wealth + moscow + poldiscST,
                     family = "binomial", data = pooled)

# models with alternative activism variable and alternative Putin support
mH3_robset8_1 <- glm(subset4 ~ openST + conscST + extraST + agreeST + neurST,
                     family = "binomial", data = pooled)

mH3_robset8_2 <- glm(subset4 ~ openST + conscST + extraST + agreeST + neurST +
                       age + edu + female + wealth + moscow,
                     family = "binomial", data = pooled)

mH3_robset8_3 <- glm(subset4 ~ openST + conscST + extraST + agreeST + neurST +
                       age + edu + female + wealth + moscow + poldiscST,
                     family = "binomial", data = pooled)


tab_model(mH3_robset5_1, mH3_robset5_2, mH3_robset5_3, 
          mH3_robset6_1, mH3_robset6_2, mH3_robset6_3,
          mH3_robset7_1, mH3_robset7_2, mH3_robset7_3,
          mH3_robset8_1, mH3_robset8_2, mH3_robset8_3,
          show.se = TRUE, show.ci = FALSE, 
          p.style = "stars", collapse.se = TRUE, show.ci50 = TRUE, show.aic = TRUE,
          file = "Table_H6.html")


### Table H7: multinomial models on 2013 data

data1$subset_mult <- ifelse(data1$putin_sup==0 & data1$polact1==1,"op_act",
                                    ifelse(data1$putin_sup==1 & data1$polact1==0,"pro_nonact",
                                           ifelse(data1$putin_sup==0 & data1$polact1==0,"op_nonact",
                                                  ifelse(data1$putin_sup==1 & data1$polact1==1,"pro_act",NA))))

data1$subset_mult <- relevel(factor(data1$subset_mult), ref = "pro_nonact")


mH3_robset3_1 <- nnet::multinom(subset_mult ~ openST + conscST + extraST + agreeST + neurST,
                              data = data1)
tab_model(mH3_robset3_1, show.se = TRUE, show.ci = FALSE, 
          p.style = "stars", collapse.se = TRUE, show.ci50 = TRUE, show.aic = TRUE,
          file = "Table_H7a.html")

mH3_robset3_2 <- nnet::multinom(subset_mult ~ openST + conscST + extraST + agreeST + neurST +
                                age + edu + female + wealth + moscow,
                              data = data1)
tab_model(mH3_robset3_2, show.se = TRUE, show.ci = FALSE, 
          p.style = "stars", collapse.se = TRUE, show.ci50 = TRUE, show.aic = TRUE,
          file = "Table_H7b.html")

mH3_robset3_3 <- nnet::multinom(subset_mult ~ openST + conscST + extraST + agreeST + neurST +
                                age + edu + female + wealth + moscow + poldisc + econ + private,
                              data = data1)

tab_model(mH3_robset3_3, show.se = TRUE, show.ci = FALSE, 
          p.style = "stars", collapse.se = TRUE, show.ci50 = TRUE, show.aic = TRUE,
          file = "Table_H7c.html")

##### Supplementary Appendix J: Redefine NA values  #####
### Missing values on leader support in Hypothesis 3 ###

# redefine NA values
data1$putin_sup_alt <- ifelse(data1$putin_sup_raw==1,1,
                              ifelse(data1$putin_sup_raw==97 | data1$putin_sup_raw==99,0,NA))
data1$putin_sup_alt2 <- ifelse(data1$putin_sup_raw==1,1,0)

# create subset variable within 2013 survey
data1$subset_alt <- ifelse(data1$putin_sup_alt==0 & data1$polact1==1,1,
                           ifelse(data1$putin_sup_alt==1 & data1$polact1==0,0,NA))
data1$subset_alt2 <- ifelse(data1$putin_sup_alt2==0 & data1$polact1==1,1,
                            ifelse(data1$putin_sup_alt2==1 & data1$polact1==0,0,NA))

### Table J1

# with first alternative (only NAs)
mH3_robset10_1 <- glm(subset_alt ~ openST + conscST + extraST + agreeST + neurST,
                      family = "binomial", data = data1)

mH3_robset10_2 <- glm(subset_alt ~ openST + conscST + extraST + agreeST + neurST +
                        age + edu + female + wealth + moscow,
                      family = "binomial", data = data1)

mH3_robset10_3 <- glm(subset_alt ~ openST + conscST + extraST + agreeST + neurST +
                        age + edu + female + wealth + moscow + poldisc + econ + private,
                      family = "binomial", data = data1)

tab_model(mH3_robset10_1, mH3_robset10_2, mH3_robset10_3,
          show.se = TRUE, show.ci = FALSE, 
          p.style = "stars", collapse.se = TRUE, show.ci50 = TRUE, show.aic = TRUE,
          file = "Table_J1.html")

### Table J2

# with second alternative (NAs added to oppositionists)
mH3_robset10_4 <- glm(subset_alt2 ~ openST + conscST + extraST + agreeST + neurST,
                      family = "binomial", data = data1)

mH3_robset10_5 <- glm(subset_alt2 ~ openST + conscST + extraST + agreeST + neurST +
                        age + edu + female + wealth + moscow,
                      family = "binomial", data = data1)

mH3_robset10_6 <- glm(subset_alt2 ~ openST + conscST + extraST + agreeST + neurST +
                        age + edu + female + wealth + moscow + poldisc + econ + private,
                      family = "binomial", data = data1)

mH3_robset10_7 <- glm(subset_alt2 ~ openST + conscST + extraST + agreeST + neurST +
                        age + edu + female + wealth + moscow + extraST*agreeST,
                      family = "binomial", data = data1)

tab_model(mH3_robset10_4, mH3_robset10_5, mH3_robset10_6, mH3_robset10_7,
          show.se = TRUE, show.ci = FALSE, 
          p.style = "stars", collapse.se = TRUE, show.ci50 = TRUE, show.aic = TRUE,
          file = "Table_J2.html")


##### Supplementary Appendix K: Robustness checks for Interaction effects #####

### Figure K1

# prepare marginal effects tables (requires output from script "Main analyses")
me.H4a <- ggeffects::ggeffect(mH3_set1_2, terms = c("extraST", "agreeST"))
me.H4b <- ggeffects::ggeffect(mH3_set1_4, terms = c("extraST", "agreeST"))

# merge marginal effects tables
me.H4a$group2 <- as.factor(rep.int(c("-1 sd", "mean", "+1 sd"), 12))
me.H4a$sample <- "2013 survey"

me.H4b$group2 <- as.factor(rep.int(c("-1 sd", "mean", "+1 sd"), 13))
me.H4b$sample <- "pooled: 2013/Navalny"

me.H4 <- rbind(me.H4a, me.H4b)

# plot
me.H4$group2 <- factor(me.H4$group2, levels = c("-1 sd", "mean", "+1 sd"))
me.H4$sample <- factor(me.H4$sample, 
                       levels = c("2013 survey",
                                  "pooled: 2013/Navalny"))

ggplot(me.H4[me.H4$group2!="mean",]) +
  geom_smooth(aes(x=x, y=predicted, color = group2),
              size=.5, se = F) +
  geom_ribbon(aes(x=x, ymin = conf.low, ymax = conf.high, fill = group2), 
              alpha = .2) +
  scale_color_discrete(name="agreeableness") +
  scale_fill_discrete(name="agreeableness") +
  xlab("extraversion") +
  ylab("predicted probability of being\nin opp. activist group") +
  facet_wrap(~sample, 
             ncol = 3) +
  theme_bw()

ggsave("Figure_K1.pdf", width = 7, height = 4) 


##### Supplementary Appendix N: Extending beyond Russia  #####

### Table N1

# m3: opposition protest (case-control design)
mH1_robset3_1 <- glm(actopp ~ openST + conscST + extraST + agreeST + stableST,
            family = "binomial", weights = weight,
            data = BYdata)

mH1_robset3_2 <- glm(actopp ~ openST + conscST + extraST + agreeST + stableST
            + age + edu + female + wealth + minsk,
            family = "binomial", weights = weight,
            data = BYdata)

mH1_robset3_3 <- glm(actopp ~ openST + conscST + extraST + agreeST + stableST
            + age + edu + female + wealth + minsk
            + private + econlastyear,
            family = "binomial", weights = weight,
            data = BYdata)

tab_model(mH1_robset3_1, mH1_robset3_2, mH1_robset3_3, 
          show.se = TRUE, show.ci = FALSE, 
          p.style = "stars", collapse.se = TRUE, show.ci50 = TRUE, show.aic = TRUE,
          file = "Table_N1.html")

### Table N2

# m1: previous protest
mH1_robset4_1 <- glm(earlyprot ~ openST + conscST + extraST + agreeST + stableST,
            family = "binomial", weights = weight,
            data = BYdata)

mH1_robset4_2 <- glm(earlyprot ~ openST + conscST + extraST + agreeST + stableST
            + age + edu + female + wealth + minsk, 
            family = "binomial", weights = weight,
            data = BYdata)

mH1_robset4_3 <- glm(earlyprot ~ openST + conscST + extraST + agreeST + stableST
            + age + edu + female + wealth + minsk 
            + luka + polint, 
            family = "binomial", weights = weight,
            data = BYdata)

tab_model(mH1_robset4_1, mH1_robset4_2, mH1_robset4_3, 
          show.se = TRUE, show.ci = FALSE, 
          p.style = "stars", collapse.se = TRUE, show.ci50 = TRUE, show.aic = TRUE,
          file = "Table_N2.html")

### Table N3

mH1_robset5_1 <- glm(luka ~ openST + conscST + extraST + agreeST + stableST,
            weights = weight,
            family = "binomial",
            data = BYdata)

mH1_robset5_2 <- glm(luka ~ openST + conscST + extraST + agreeST + stableST
            + age + + edu + female + wealth + minsk,
            weights = weight,
            family = "binomial",
            data = BYdata)

mH1_robset5_3 <- glm(luka ~ openST + conscST + extraST + agreeST + stableST
            + age + + edu + female + wealth + minsk
            + private + econlastyear,
            weights = weight,
            family = "binomial",
            data = BYdata)

tab_model(mH1_robset5_1, mH1_robset5_2, mH1_robset5_3, 
          show.se = TRUE, show.ci = FALSE, 
          p.style = "stars", collapse.se = TRUE, show.ci50 = TRUE, show.aic = TRUE,
          file = "Table_N3.html")

### Figure N1

# predict
predictBY_H3_ex <- ggeffects::ggeffect(mH1_robset3_2, terms = c("extraST [meansd]"))
predictBY_H3_ag <- ggeffects::ggeffect(mH1_robset3_2, terms = c("agreeST [meansd]"))

# plot
pBY_H3_extra <- ggplot(predictBY_H3_ex[predictBY_H3_ex$x!=0.16,]) +
  geom_pointrange(aes(x = factor(x), y = predicted,
                      ymin = conf.low, ymax = conf.high)) +
  geom_text(aes(x = factor(x), y = predicted, label = round(predicted, 2)),
            size=4, hjust = 1.3) +
  scale_x_discrete(labels = c("-1 sd", "+1 sd")) +
  xlab("extraversion") +
  ylab("predicted probability of being\nin opp. activist group") +
  ylim(0.18, 0.8) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)); pBY_H3_extra

pBY_H3_agree <- ggplot(predictBY_H3_ag[predictBY_H3_ag$x!=0.06,]) +
  geom_pointrange(aes(x = factor(x), y = predicted,
                      ymin = conf.low, ymax = conf.high)) +
  geom_text(aes(x = factor(x), y = predicted, label = round(predicted, 2)),
            size=4, hjust = 1.3) +
  scale_x_discrete(labels = c("-1 sd", "+1 sd")) +
  xlab("agreeableness") +
  ylab("\n") +
  #ggtitle("agreeableness") +
  ylim(0.18, 0.8) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)); pBY_H3_agree

ggsave("Figure_N1.pdf", gridExtra::grid.arrange(pBY_H3_extra, pBY_H3_agree, 
                                            ncol = 2, top = "Belarus survey, case-control design"),
       width = 7, height = 3.5)

